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ABSTRACT 

We study oscillations of slowly rotating relativistic barotropic as well as non-barotropic 
polytropic stars in the Cowling approximation, including first order rotational correc- 
tions. By taking into account the coupling between the polar and axial equations, we 
find that, in contrast to previous results, the m = 2 r modes are essentially unaf- 
fected by the continuous spectrum and exist even for very relativistic stellar models. 
We perform our calculations both in the time and frequency domain. In order to nu- 
merically solve the infinite system of coupled equations, we truncate it at some value 
'max- Although the time dependent equations can be numerically evolved without any 
problems, the eigenvalue equations possess a singular structure, which is related to the 
existence of a continuous spectrum. This prevents the numerical computation of an 
eigenmode if its eigenfrequency falls inside the continuous spectrum. The properties 
of the latter depend strongly on the cut-off value Z ma x and it can consist of several 
either disconnected or overlapping patches, which are the broader the more relativistic 
the stellar model is. By discussing the dependence of the continuous spectrum as a 
function of both the cut-off value l max and the compactness M/R, we demonstrate 
how it affects the inertial modes. Through the time evolutions we are able to show 
that some of the inertial modes can actually exist inside the continuous spectrum, 
but some cannot. For more compact and therefore more relativistic stellar models, the 
width of the continuous spectrum strongly increases and as a consequence, some of 
the inertial modes, which exist in less relativistic stars, disappear. 

Key words: relativity - methods: numerical - stars: neutron - stars: oscillations - 
stars: rotation 



1 INTRODUCTION 

Non-radial oscillations of relativistic stars gained a lot of interest in the last decades because of the possible detection of 
their associated gravitational waves. Especially after the discovery that the r modes of rotating neutron stars are generically 
unstable to gravitational radiation (Andersson 1998, Friedman & Morsink 1998), the interest in the rotational instabilities 
renewed. Great effort has been put into understanding this instability and its implications on the spin evolution of both newly 
born and old accreting neutron stars. If not damped by viscous effects, the r mode instability should be a very efficient way 
to generate gravitational radiation. For a recent account of the r-mode instability we refer to the review by Andersson and 
Kokkotas (2001). 

Up to now, most of the studies of the r-mode instability have been done within the Newtonian framework. In a first 
attempt to study the relativistic analogue, Kojima (1997, 1998) considered the axial equations for slowly rotating stars in the 
so-called "low-frequency approximation". He found a situation quite different from Newtonian theory, in that to first order 
the relativistic frame dragging in rotating stars gives rise to a continuous spectrum of frequencies instead of a single r-mode 
frequency. In a subsequent work, Beyer & Kokkotas (1999) proved in a mathematically rigorous way the existence of this 
continuous spectrum in the low-frequency approximation. 

Lockitch, Andersson & Friedman (2001) showed that for uniform density stars, in addition to the continuous spectrum, 
Kojima's equations also admits discrete mode solutions, the relativistic r modes. Ruoff & Kokkotas (2001) and Yoshida (2001) 
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extended this study to polytropic and realistic models and showed that for certain stellar models, no r modes can be found, 
even when the low-frequency assumption is dropped and the radiation reaction is included (Ruoff & Kokkotas 2002; Yoshida 
& Futamase 2001). In all these studies, however, only the axial perturbations have been considered, and any coupling to the 
polar equations has been neglected. Lockitch et al. (2001) have argued that this assumption is only justified for non-barotropic 
stars, as the relativistic r modes in barotropic stars should be hybrids, i.e. they retain a non-vanishing polar contribution in 
the non-rotating limit. 

By including higher order rotational corrections and the coupling to the polar equations, Kojima & Hosonuma (2000) 
derived a fourth order equation, an extension of Kojima's master equation, which as a consequence does not have the singular 
structure of the latter. Although not rigorously proved, it was argued that this could possibly resolve the problem of the 
non-existence of the r modes in certain stellar models (Lockitch & Andersson 2001). 

Recently Yoshida & Lee (2002) announced that they can find relativistic r modes in the Cowling approximation for 
non-barotropic stellar models. They obtained their results by only including first order rotational correction terms. However, 
they did not restrict themselves to merely considering the coupling between the axial equation of order I to the polar equations 
of order I + 1, but instead they solved the complete system starting from / = jm| up to some value i m ax, which they chose 
sufficiently high enough to ensure the convergence of the modes. Lockitch et al. (2001) had already formulated this problem 
in a somewhat different way, by neglecting only the perturbations which vanish in the stationary limit, such as the energy 
and pressure perturbations for instance, and some of the metric variables. So far they have presented numerical results only 
in a post Newtonian treatment of relativistic r modes for uniform density stars. Their numerical eigenfunctions were used by 
Stergioulas & Font (2001) as approximate r-mode initial data for a 3-D general-relativistic hydrodynamical time evolution in 
the Cowling approximation. They found that the r mode oscillates at a predominately discrete frequency with an indication 
of the kinematical drift, as predicted by Rezzolla, Lamb & Shapiro (2000). 

In this paper, we follow the approach of Yoshida & Lee (2002) and extend it also to the time domain. Furthermore, we 
consider both barotropic and non-barotropic stellar models and discuss in detail how the continuous spectrum appears and 
how it is depends on the cut-off Z max . By computing a sequence of stellar models with different M/R ratios, we study how 
the continuous spectrum affects their inertial modes. In this paper, we call all the modes "inertial modes", whose restoring 
force is mainly the Coriolis force. Thus we consider the r modes as a subclass of inertial modes whose motion is dominated 
by the axial velocity field. 

The paper is organized as follows: In Section 2 we present both the time dependent and the eigenvalue forms of the coupled 
perturbation equations. We first specialize on the axisymmetric case, where we describe how to compute the continuous 
spectrum and show its dependence on / max . The numerical methods and results are presented in Section 3 and we conclude 
with the discussion in Section 4. Throughout the paper we use geometrical units G = c = 1. The over-dots and primes denote 
differentiation with respect to the time and the radial coordinates, respectively. All perturbation variables should carry indices 
I and m, however, we conveniently omit them. Only in cases, where the I index is important, we will include it. Finally, we 
assume m > 0, as one can always transform a m < mode with frequency a into a m > mode with frequency —a by taking 
its complex conjugate. 



2 FORMULATION OF THE PROBLEM 
2.1 Background model 

We consider a uniformly rotating relativistic star with angular velocity Q, which is assumed to sufficiently small that any 
deviation from sphericity can be neglected. The metric can then be written as (Hartle 1967) 

ds 2 = ~e v dt 2 + e 2X dr 2 + r 2 {d6 2 + sin 2 6d(j> 2 ) - 2lo sin 2 6dtd(j>. (1) 
Assuming a perfect fluid star, the energy momentum tensor 

T^ u = (p + e) u^Uv + pg^, (2) 
with pressure p, energy density e and four-velocity components 

u — e~", (3) 
it* = ne~". (4) 
We will use a relativistic polytropic equation of state given by 

P = *p r , (5) 
e = P+f^T- (6) 

where k denotes the polytropic constant, F is the polytropic index and p the rest mass density. With this form of the equation 
of state, the polytropic index F is equivalent to the adiabatic index and obeys the relation 
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r _ p + e dp 
p de 



2.2 Perturbed equations of motion 

We assume the oscillations to be adiabatic, so that the relation between the Eulerian pressure perturbation Sp and energy 
density perturbation Se is given by 

p-+-e 6e+p t It- 1 )' 

where Fi represents the adiabatic index of the perturbed configuration and £ r is the radial component of the fluid displacement 
vector £ M . The sound speed C s is 

2 Ti dp 



dp 



(8) 



C: 



T de' 



0) 



The complete set of the perturbed Einstein equations has been derived using the BCL gauge (Battiston, Cazzola & Lucaroni 
1971) in Ruoff, Stavridis & Kokkotas (2001), which we will refer to as RSK. As we are working in the Cowling approximation, 
we only need to consider the five fluid perturbations, which are the three components of the velocity perturbations Sui, the 
enthalpy perturbation H and the radial component of the displacement vector £ M . Following RSK, we expand these quantities 
as 



Sun — 



5e 



e = 



-e >«i Yim, 

1 ,777 
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l, in 
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(p + e) 2 

Tip 
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(10) 

(ii) 

(12) 
(13) 
(14) 



Later we will introduce new variables «2 and U3 instead of u 2 and H3. These include rotational corrections and therefore differ 
from the definitions of u 2 and 113 in RSK. The time dependent perturbation equations follow directly from Eqs. (68)-(72) of 
RSK with all the metric perturbations set to zero: 



(8t + irnQ)H = (';<),'" 
(dt + imQ) tti 
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1 ip 
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(d t + imQ) 112 — H + (imu 2 + C^us) — A, 

(dt + imQ) u :i = 

(dt + imQ,) £ = 
where 

A = voC 2 e~ 2X u 1 + I 'lu — A H — I ui — e~ '— -U2 + w \ v 1 + lu e ~"ui 

\ rj r 1 J L V rj 

and 

A = 1(1 + 1). 

The operators -C* 1 , C^ 1 , and /Z^ 1 are the same as in RSK and are defined by their action on a perturbation variable P 
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(15) 

(16) 

(17) 
(18) 
(19) 

(20) 



£-±\ pirn. 



(I - l)Qi m P l ' lm -(1 + 2)Q l+lm P 



l + lm 



C^P lm = -(l + l)Qi m P'- lm + lQi + i m P 

£±1 plm ( j i\/f 1 1 \ n j-jl — lm 



l + lm 
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with 

Qlm 



(I -m)(l + m) 



(25) 



(21 - 1)(2Z + 1) ' 

By changing to new variables, we can simplify the above set of evolution equations. To this end we use Eq. ( |l5| ) and rewrite 
Eq. IH as 

A = voe~ 2v d t H + e~ 2X Bm + 0(Q. 2 ) (26) 
with 

r 



B 



Id + 2vj ( V 



(27) 



Neglecting second order terms and moving the terms including d t H to the lefthand side, we can write Eqs. |t7| ) and ( |l8| ) as 

(28) 
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This suggests to define new variables 



U2 := U2 + vm — e fi, 



t*3 := 



(29) 



(30) 
(31) 



and rewrite Eqs. (|l^)-(p^|) in terms of 112 and 113. After discarding any second order terms introduced by this replacement the 
perturbations equations now read 
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(32) 
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(34) 
(35) 
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In order to find the eigenmodes of this coupled system of equations, we assume our perturbations variables to have a harmonic 
time dependence exp(iat). By replacing all time derivatives by icr and letting H — > iff, 113 — > k/,3 and £ — > i£ in order to obtain 
a purely real valued system of equations, we obtain two ODEs for H and mi and three algebraic relations for m 2 , M3 and £. 
The relation for £, which follows from Eq. (E(3) is rather simple and can be immediately used to eliminate £, giving 
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where we have defined 



E := a + mfi — 
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(37) 
(38) 
(39) 
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The relations (Q and © cannot be solved directly for U2 and M3 as they involve the coupling operators Cf 1 . We should 
stress again that because of this coupling, both sets of equations (^)-(^) and j3^-d4C,) have to be considered as infinite 
systems with I running from m to infinity. There is no coupling to equations with I < m, since from its definition ([^ it 
follows that the coupling coefficient Q mm , which would couple to I = m — 1, is zero. Furthermore, the equations form two 

© 2002 RAS, MNRAS 000, §4l| 



Inertial modes of slowly rotating relativistic stars in the Cowling approximation 5 



independent sets, each belonging to a different parity. In the non-rotating case, one usually distinguishes between polar and 
axial perturbations. Under parity transformation polar perturbations change sign as ( — 1)', axial perturbations as (— In 
the rotating case, however, the polar and axial equations are coupled, and therefore the distinction between polar and axial 
modes cannot be upheld any more. However, the equations do not mix the overall parity as can be seen as follows. A polar 
equation with even / = m has even parity. It is coupled to an axial equation with I + 1, whose parity is also even. The next 
coupling is again to a polar equation with I + 2, thus having even parity. As this continues in the same manner, this means 
that for even m, the complete coupled system with a leading polar equation has even parity. Conversely, the other system 
starting with a leading axial equation has odd parity. For odd m, we obtain the reversed situation. This implies that for any 
given m, it makes sense to distinguish the modes according to their overall parity. Lockitch & Friedman (2000) introduced 
the notion of polar or axial led modes, depending on whether the leading equations with / = m are axial or polar. For even 
m, the polar led modes have even parity and the axial modes odd parity, and vice versa for odd m. The axisymmetric case 
m = is somewhat special in that only the even parity modes start with I = whereas the odd parity modes have to start 
with I = 1, as there are no axial I = equations. In the following section, we will restrict ourselves to positive m. 



2.3 The Continuous Spectrum 

For simplicity, we now focus on axisymmetric perturbations with m — 0, which leads to considerable simplifications of the 
time independent equations ©-©. With the expansion of the operators C^ 1 according their definitions ( p2j ) - the 
m = equations read (we now explicitly include the index I): 
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For / = 0, we have = u% = 0. To solve the ODEs (H) and (j43|), we have to compute it| and u l 3 in terms of H l and u[. 
This can be done the most easily by combining u l 2 and u 3 as well as H l and u[ into vectors and rewriting the two algebraic 
relations (^) and |f5| ) as 
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Defining 
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we can cast Eq. (^?j) into a more compact form 
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This, again, can be viewed as a matrix equation for two infinitely dimensional matrices U and S acting on the vectors 
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whose respective elements are the 2- vectors u l and s l . Both U and S are tridiagonal block matrices, with each block given by 
the above 2x2 matrices. Explicitly, U is given by 

/ U„ U\ ... \ 

ul u a ul 
ul 

(54) 
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A similar structure holds for S. Now we can write Eq. (M) as 



which can be solved for u l by multiplying both sides with U 1 

1 = 1... oo. 
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It is important to note that the matrix U is r-dependent since its elements U± contain the function w. This raises the question 
whether one can assure that U will be invertible for any frequency a at each value of r in the stellar interior. For U to be 
invertible its determinant must not vanish. 

In the Newtonian limit w — ► Q, the matrix U does not depend on r any more, and one can easily show that U becomes 
singular only for a discrete set of values of a. These frequencies represent the solutions of the homogeneous part Um = of 
Eq. (pB|). In the relativistic case, however, we have to replace Q. by the r-dependent effective angular velocity zu. Thus the 
zeroes of det U now depend on the position inside the star. This means that each single Newtonian frequency will be spread 
out in the relativistic case into a continuous band of frequencies, determined by the values of zu at the centre and at the 
surface of the star, which we will denote by zuo and zur, respectively. The total range of the continuous spectrum is then the 
sum of all the individual bands originating from their Newtonian discrete values. 

In order to numerically perform the inversion and solve equation ([i(]), we have to truncate the system at some value 
( max > 1. Let us consider the simple case, where we truncate the system at Z max = 2, i.e. we couple only / = 1 with / = 2 and 
neglect all I > 3. Then the matrix U is given by 



U a 

Ul 



Ua 



a 








—dzu 





(7 


— &ZUQ2 








-2roC 


h <r 





2zuQ 2 








a 



Q2 \ 



with Q2 = 1/v 15 and det U = yields 

2 4 2 n o 2 

a — — zu — U.o zu . 

5 

Hence, the range of the continuous spectrum is given by 
0.894 ceo < cr c < 0.894 zu R , 

if we restrict ourselves to positive frequencies. If we also include the I = 3 terms, we obtain 

2 84 2 ^ , „, , 2 

a = — zu w 1.714 zu . 
49 

and the range is shifted to 
1.309 too < a c < 1.309 wr. 

Going one step further and including / = 4, we obtain two solutions 



4 
441 



(147 ±42^7) vj 2 , 
corresponding to the frequency intervals 
1.530 zuq < <7 C < 1.530 zur. 
for the plus sign and 
0.5704to < <t c < 0.5704 zu R 



(57) 

(58) 
(59) 
(60) 
(61) 
(62) 
(63) 
(64) 
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for the minus sign. Adding I = 5 will still yield two ranges and for I — 6 we obtain three. The pattern is as follows: When a 
coupling to a higher even value of I is included, one additional range of the continuous spectrum will appear. Depending on 
the respective values of vao and vjr, the various ranges might as well overlap. This actually happens for the more relativistic 
stars since the variation in w is much greater. 

A similar picture also holds in the case Here, each coupling to a higher value of / results in one additional band of 

the continuous spectrum. We have already stated that in the case of the matrix U becoming singular at some point inside the 
star, we cannot invert it anymore and the integration of the eigenvalue system (p^)-(ffO|) fails. In the purely axial case, Ruoff 
& Kokkotas (2001) showed that by a series expansion one could still obtain mathematically valid mode solutions inside the 
continuous spectrum, however, only at the price of a having a divergence in the associated fluid perturbation at the singular 
point. Thus they concluded that these modes are unphysical and should be discarded. The question is whether we can adopt 
this point of view and claim that inside the continuous bands modes cannot exist. 

It is clear that the widths of the continuous bands strongly depend on the stellar parameters and go to zero in the 
Newtonian limit. For weakly relativistic stellar models, each individual band is fairly small, but the total range, which is 
the sum of these small bands, can become quite large. If it is true that modes should not exist inside the continuous bands, 
then even for weakly relativistic stellar models, there should be almost no true mode solutions. As we cannot perform the 
integration of ©-(0) inside the continuous bands, we cannot tell whether there still might be valid mode solutions. Trying 
series expansion of the coupled equations around the singular point will be quite involved. This is why we rather rely on the 
evolution of the time dependent equations, which are free of singularities. By examining the Fourier spectra we can assess 
whether or not we can find modes inside the continuous spectrum. As we shall show and discuss, there are indeed modes 
which lie inside the continuous bands. However, as the stellar models become more compact, the continuous spectrum has 
the ability to destroy more and more of the inertial modes. 



3 NUMERICAL METHODS AND RESULTS 

The numerical evolution of the coupled equations does not present any severe numerical problems. Although one might 
encounter long term instabilities, if a large number of I is included, these can usually be overcome by increasing the spatial 
resolution or shifting the numerical origin some grid points away from zero. Up to values of about Z max = 10, evolutions with 
200 grid points inside the star yield quite accurate results for the inertial modes. As the inertial modes are located at at 
the lower end of the frequency spectrum, we have to perform long time evolutions in order to obtain a sufficient frequency 
resolution. Usually, we have to evolve up to the order of 500 ms, corresponding to a total number of time steps well above 10 6 . 

The numerical mode calculation is somewhat more troublesome. As explained in the previous section, we cannot compute 
the modes once they are inside the continuous spectrum, because the matrix U becomes singular and inversion is no longer 
possible. To find the modes outside the continuous spectrum, we proceed as follows. We first rescale the variables according 
to 

H l -> r l H l , (65) 
u-l —> r it 1; (do) 
U2 — » r«2, (67) 
u l 3 -» r l+1 u 3 . (68) 

This ensures that all our new variables have Taylor expansions around the origin starting with the constant term. To initiate 
the integration from the centre toward the stellar surface, we have to prescribe initial data for each H l and u\. At each 
integration step, we have to invert the matrix U in order to compute u\ and M3. At the origin, H l and u\ are not independent, 
but are related to leading order by 

H (0) = \a + rntt - 4uj i: l __ 1 ^ Q lm I — - — . (69) 

To leading order, however, there is no relation between H l and any other H l , so that we have n — / max — m + 1 degrees of 
freedom for the initial values, as I runs from m to Z max . The frequency a represents a further degree of freedom, which we 
have to choose deliberately. The missing boundary conditions, which fix these degrees of freedom, come from the requirement 
of a vanishing Lagrangian pressure perturbation at the stellar surface. Since this requirement translates into the simultaneous 
vanishing of the quantities 



Ap l 



(a + rml)e n — v U\ 



, m < I < /max, (70) 

we obtain Z max — m + 1 conditions. We have thus reduced the number of degrees of freedom to one, which is the value one 
expects for linear systems, as the overall amplitude of the perturbations remains arbitrary. Let us denote a set of initial data 
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for the H l by Tik := \H k n , . . . , H k '" ax ^ . In order to compute the modes, we choose a frequency a and perform n integrations for 
the following n different sets of initial data: Tri\ = (1, 0, . . . , 0) , TL2 = (0, 1, 0, . . . , 0) , . . . , TL n = (0, . . . , 1). The respective initial 
data for ui follow from (|69[). For each set TCk, we compute the n Langragian pressure perturbations Ap l k , which in general will 
not be zero. After having performed the n integrations, we end up with n x n values for the Langrangian pressure perturbations 
Ap l k , which can be arranged inanxn "pressure matrix" P. As the n solutions pertaining to the n sets of initial data TLk are 
linearly independent, we can try to use them to construct another solution by a suitable linear combination, which satisfies 
Ap l = for all I simultaneously. In other words we have to find n coefficients ait, such that the linear combinations a k Ap l k 
vanish for all I. This is possible if and only if the determinant of the pressure matrix P vanishes. The algorithm for finding 
the eigenmode frequency is then obvious. For each a, we construct the pressure matrix P and compute its determinante. If 
zero, we have found an eigenmode. 

To also compute the eigenfunction, we have to solve the homogeneous system of equations 

'max 

^ a k Ap k =0, l=m... /max (71) 

k — m 

for the coefficients a k , which can be easily done by a "singular value decomposition" of P (see e.g. Press et al. 1992). The 
eigenfunction can be obtained by the linear combination of the n solutions corresponding to the n sets of initial data Tik, 
with the appropriate weighing coefficients ak- 

Although this method is quite easy to implement, it suffers from the drawback of breaking down when i max > m + 7. 
This failure results from the insufficient prescription of the initial data, where we should include higher order terms in the 
Taylor expansion, as the leading orders are set to zero for all I but one. For the purpose of computing the inertial modes for 
small values of I, however, our method is robust enough, and we obtain perfect agreement with the frequencies obtained from 
the time evolutions. Different methods are used by Lockitch et al. (2001) and Yoshida & Lee (2002). 

We have checked the consistency of our codes by feeding the evolution code with the eigenfunctions obtained from the 
eigenvalue code. As expected the time evolution yields a standing wave with one single frequency. Arbitrary initial data excite 
the / and p modes together with the inertial and in the non-barotropic case also the g modes. The frequencies found by 
Fourier transforming the time sequences agree with those found from the mode calculation within less than one per cent. This 
difference depends only on the numerical resolution of the time evolutions and tends to zero for increasing resolution. 



3.1 The axisymmetric case 

Let us first assess the influence of the rotation on the / and p modes. We take the same stellar model as Font et al. (2001), 
which is a T = 2 polytrope with k = 217.86 km 2 . (Note their value k = 100 is in dimensionless units where G = c = Mg = 1). 
With a central rest mass density po = 7.914 x 10 14 g/cm 3 , the total mass is M = 1.4 Mq and the radius is R = 14.15 km. 
The rotation rate is given by Q = 4230 s -1 corresponding to e — 0,/Q.k = 0.5223, where Qk = \J M/R 3 . In the following 
we consider only barotropic perturbations, i.e. we use Ti = 2, thus suppressing the g modes. We will study non-barotropic 
perturbations in the non-axisymmetric case. 

In Table 1, we show the change in frequencies of the (quasi) radial f° and p° modes, when the coupling to higher values of 
I is included. As we are considering axisymmetric perturbations, there is no m-splitting and rotational effects only come from 
the coupling of the equations with different I. For Z max = 0, there are no rotational corrections, i.e. the equations describe the 
non-rotating case. For Z max = 1, the polar I = equations couple to the axial I = 1 equation for 113, and we find a shift in the 
frequencies towards higher values. This shift is about 3.7 per cent for the f° mode, 1.3 per cent for the p\ mode and even less 
for higher modes. When further couplings to higher I are included, the additional frequency corrections become smaller and 
the frequencies converge rapidly. Beyond i max = 3, the frequency corrections are less than 0.01 per cent. Hence, we see that 
the quasi-radial / and p modes are quite insensitive to rotation. The effects of rotation on the axisymmetric modes have been 
investigated using the fully nonlinear equations in the Cowling approximation by Font et al. (2001). However, a comparison 
between our results and theirs is not very meaningful as they study rapidly rotating models, which undergo quite large radius 
and mass corrections with respect to the non-rotating case. In the perturbation formalism, these corrections appear only at 
the Q, 2 level, so that we cannot account for them when only first order corrections are included. That these corrections are 
important can be seen from the fact that their frequencies decrease with increasing rotation rate while ours actually increase. 
By including some of the second order corrections we have convinced ourselves that the frequencies get pushed down to smaller 
values, which then agree much better with the ones of Font et al. (2001). Therefore, we compare our values in Table 1 with 
theirs only in the non-rotating limit. For the first four radial modes, they find the frequencies 2708 Hz, 4547 Hz, 6320 Hz and 
8153 Hz, which agree very well with the first row of Table 1. For m ^ 0, we find the well-known linear splitting of the / and 
p modes into a lower and a higher frequency part. 

Let us now turn our attention towards the inertial modes. In the axisymmetric case, the first inertial modes appear, when 
coupling to I — 2 is included. As described in the previous section, a continuous band of frequencies appears at the same time. 
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Table 1. f° and p° mode frequencies as functions of Z max for the 1.4Mq model described in the text. The rotation rate is e = 0.5223. 
The frequencies are given in Hz. The / and lowest p mode experience the strongest rotational corrections. For a coupling up to i max = 3, 
the mode frequencies have converged. 



^max 


f° 


P? 


Pi 


Pi 


Pi 





2687 


4551 


6344 


8111 


9867 


i 


2787 


4610 


6386 


8144 


9894 


2 


2795 


4613 


6387 


8145 


9895 


3 


2796 


4613 


6387 


8145 


9895 



In the time evolution, the modes can be distinguished from the continuous band as follows. Place some observers at different 
locations inside the star. Take the FFT of each observer's time series and plot the various spectra on top of each other. Those 
peaks that coincide for all observers belong to modes, while those peaks which are different for each observer correspond to 
frequencies from the continuous spectrum. 

As discussed above, the location of the continuous spectrum depends strongly on how many Z are taken into account. 
Depending on Z max , it can consist of several patches. The width of each patch is determined by the compactness of the stellar 
models and shrinks to a set of discrete values in the Newtonian limit. For very compact stars, the patches can become quite 
broad and in general overlap. In the following, we will investigate the behaviour of both the continuous spectrum and the 
inertial modes as functions of the stellar compactness M/R and the value of Z max . As a first step, however, it is helpful to 
understand the behaviour of the inertial modes for an almost Newtonian stellar model, where the ranges of the continuous 
spectrum are very narrow. 

Lindblom & Ipser (1999) have shown that for the Newtonian Maclaurin spheroids, the eigenvalue problem can be made 
separable by choosing appropriate spheroidal coordinates. Hence, Z and m are "good" quantum numbers, which can be used to 
label the modes. For each value of m > and I > m there exist l — m inertial modes. For m — 0, one finds Z — 1 mode solutions. 
In this case, however, the eigenvalue problem can be written in terms of a 2 , hence both a and —a are valid eigenfrequencies. 
For even I, there is an odd number of eigenvalues, hence one eigenvalue has to be a — (c.f. Table 1 of Lindblom & Ipser 
1999). If one restricts oneself to positive frequencies, there is no mode for I = 2, one mode for Z = 3 and 4, two modes for 
I — 5 and 6 and so on. 

In our case, I is not a good quantum number any more, but we still can assign a certain value Z to each mode as we shall 
demonstrate. To this end we now study for an almost Newtonian stellar model with compactness M/R = 0.01 how both the 
continuous spectrum and the inertial modes are affected as Z max is successively incremented. 

Figure 1 shows the locations of the continuous bands (grey shaded areas) together with the rotational modes one can find 
for different values of Z max , ranging from Z max = 2 to 5. For Z max = 2, we have a single continuous band in whose vicinity many 
modes can be found. These become more densely spaced as they approach the continuous spectrum. The modes to the left 
of the continuous spectrum have even parity and the modes to the right have odd parity. When we go to Z max = 3, we obtain 
a similar picture. The continuous spectrum has moved to the right and again is surrounded by a large number of modes. In 
addition, we find a more isolated mode closer to the position, where the continuous spectrum for Z max = 2 was located. This 
is actually the first "true" inertial mode, as it does not change if we further increment Z max . All the other modes which we can 
find for Z max = 2 and Z max = 3 are "fake" modes, i.e. they are artefacts due to the truncation of the equations and dissappear 
again when the coupling to higher Z is included. We will label the "true" inertial modes with 4?, where Zo is the value of Z max 
for which it appears for the first time and does not change its frequency significantly when Z max is further increased. Since for a 
given Z max , more than just one new inertial mode can appear, we order these modes through the lower index n with ascending 
frequency, i.e. a l ° < cr^° +1 . According to this convention, the inertial mode appearing for Z max = 3 is the if mode, which has 
odd parity. Our upper index Zo agrees with the definition of Yoshida & Lee (2000a) and differs from the one of Lockitch & 
Friedman, who use Zo — 1 to denote the same mode. Yoshida & Lee (2000a), however, use the lower index n in a somewhat 
different form, as they compute the modes in the rotating frame, which then form a set with equal fractions of negative and 
positive frequency modes. Yoshida & Lee attach negative indices n to the negative frequencies modes and positive n to the 
positive frequency modes. 

For Z max = 4, the continuous spectrum consists of two patches, and we find the i\ mode, which now has even parity. 
For Z max = 5, these two patches shift to the right and two new inertial modes, i\ and i\ appear. The emerging pattern is as 
follows: For each even Z max , a new patch of the continuous spectrum appears. If Z max > 2 is odd, (Z max — l)/2 new odd parity 
inertial modes appear, if Z max is even, we find (Z max — 2)/2 new even parity modes. Besides the inertial modes for a given Z max 
appear close to the positions where the patches of the continuous spectrum were located for Z max — 1. Examination of the 
eigenfunctions reveals that a mode i n has no node in the coefficient H , one node in the coefficient H l ~ 2 , and n nodes in the 
coefficient H l ~ 2n . 
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It can be noticed that for Z max = 4 and 5, there are only very few modes in the vicinity of the left patches of the 
continuous spectrum, which is due to the fact that these patches are somewhat broader than the other ones. If one considers 
less relativistic models, the width of the patches shrink and more modes can be found around them. In the Newtonian case, 
where all the patches degenerate to single points, we expect an infinite number of modes around these points. As the star 
becomes more relativistic and the continuous bands expand, they will swallow more and more modes. In principle this is of 
no relevance as these modes are not physical, anyway. The remaining question is, whether or not the continuous spectra are 
able to swallow the "true" inertial modes, if they grow broad enough. 

The answer is depicted in Fig. 2, where we plot the inertial modes together with the continuous spectra as functions of 
both the compactness M/R and £ max . Here, we plot only the "true" inertial modes and not those resulting from the truncation 
of the equations at Z max as in Fig. 1. The four panels illustrate the situation for the different cut-off values Z max ranging from 
3 to 6. The compactness M/R is from 0.01 (as in Fig. 1) to 0.22. The stability limit with respect to radial collapse is reached 
at M/R = 0.215. All stellar models have the same e = 0.5223, which means that the less compact models rotate much more 
slowly than the highly relativistic ones. Since the frequencies of the inertial modes depend essentially linearly on the rotation 
rate, they would become very small for less compact models. This is the reason why we express them in units of fix- We 
thus see that the inertial modes lie in the interval between and 1 for the whole compactness range. (Actually they lie in the 
interval between -1 and 1, but as negative and positive frequencies have the same magnitude, we plot only the positive ones). 

Starting with i max = 3, we have one single continuous spectrum, which broadens considerably as the compactness M/R is 
increased. For small M/R, we also find the if mode, which moves to lower frequencies as M/R increases. However, the lower 
border of the continuous spectrum decreases much faster and for M/R ~ 0.17 it reaches the mode. For larger M/R the mode 
cannot be found any more. We stress again that once the mode would be inside the continuous spectrum, the eigenvalue code 
breaks down because the matrix U becomes singular and inversion is no longer possible. In this case, we take the results of 
the time evolutions to decide whether or not the mode exists inside the continuous spectrum. And indeed, the time evolution 
shows that for M/R > 0.17 the mode does not penetrate the continuous spectrum but instead disappears. 

When we go to Z max = 4, the same happens to the if mode which reaches the upper patch of the continuous spectrum 
at about M/R ~ 0.11 and disappears for larger M/R. Surprisingly, the if mode now exists over the complete compactness 
range. Indeed, the time evolution shows that it can even exist inside the continuous spectrum for M/R > 0.2. This feature is 
even more apparent when the coupling up to i max = 5 is included. In this case, the if mode is always inside the continuous 
spectrum over the whole compactness range. However, the i\ mode still does not penetrate the continuous spectrum and the 
same is true for the if and i| modes. Notice also that for M/R > 0.17, the two patches of the continuous spectrum begin to 
overlap. 

For i max = 6, both the if and if modes exist for the complete compactness range. Again, the if is always inside the 
continuous spectrum, whereas the if mode enters it only for M/R > 0.2. The i| mode, in contrast, ceases to to exist for 
M/R > 0.1. For the / max = 7 case, which is not shown, only the if, if and if modes exist for the complete compactness range, 
most of the others exist only outside the continuous spectrum. Some actually do penetrate the continuous spectrum, but if 
the star becomes too compact, they eventually vanish. 

The emerging picture is as follows. Only the modes with low frequencies, such as the if, if and if modes for instance 
are able to exist for the complete compactness range and even inside the continuous spectrum. Those with higher frequencies 
exist only for a much smaller range. This might be attributed to the fact that the patches of the continuous spectrum are 
much broader for larger frequencies. As a consequence, the total number of inertial modes decreases as the stellar models 
become more relativistic. However, this leaves the question why some of the modes can survive inside the continuous spectrum 
and some cannot. We conjecture that this has to do with the value of Z max , where the mode occurs for the first time. If for 
a certain value of i max it appears well outside the continuous spectrum, then it will survive if i max is further increased and 
the mode suddenly finds itself inside the continuous spectrum. In other words, it is not the mode which moves inside the 
continuous spectrum as Z max is increased, it is the continuous spectrum which shifts itself to the position where the mode is 
located. If, however, for a certain Z max the continuous spectrum is already sufficiently broad to occupy the position where an 
inertial mode would appear, then this does not seem to be able to establish itself, even when higher values of / are included. 

3.2 The non-axisymmetric case 

For m 0, we find essentially the same picture as in the axisymmetric case, i.e. some of the modes can exist for the whole 
compactness range while some cannot. Here, the most interesting question, of course, is what happens to the relativistic r- 
modes and in which way are they affected by the continuous spectrum. One can see that the relativistic case differs qualitatively 
from the Newtonian one, as in the latter, one can immediately show that the leading order of the r-mode frequency (in the 
inertial frame) is given by 



a 




(72) 
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whereas the analogous relativistic calculation yields 
a(r) = -mfi 1 - / i W W 



Z(Z + 1) V fi 



(73) 



Instead of a single frequency, we have a continuous spectrum of frequencies, determined by the values of the frame dragging uj 
at the centre and surface of the star. Herein, we have chosen the sign of a such that it yields negative frequencies for positive 
values of m. 

The above values for a were obtained by assuming that the r mode is purely axial to leading order, and the coupling 
to the polar equations is just a higher order correction. However, it was shown by Lockitch et al. (2001) that this is not 
true for barotropic perturbations, as this assumption leads to some incompatibilities within the polar equations. Instead, 
the barotropic r mode necessarily has to include polar contributions, which do not vanish in the non-rotating limit. For a 
"pure" r mode, the polar contribution has to vanish. In the terminology of Lockitch et al. (2001), the relativistic barotropic 
r mode is therefore an "axial-led hybrid" mode. In the non-barotropic case, things are different, as the r mode cannot have 
non-vanishing polar contributions in the non-rotating limit. This is because for non-rotating non-barotropic stellar models, 
there are no stationary polar perturbations, since they would immediately excite g modes. For non-rotating barotropic stars 
on the other hand, the g modes form a degenerate space of zero-frequency modes, which is why the axial inertial modes can 
have polar contributions. Thus, for non-barotropic stars, the r mode should not be an axial-led hybrid, but it should reduce in 
the non-rotating limit to a purely axial stationary perturbation. Nevertheless, we will call the inertial mode i\ with I = m + 1 
an r mode, regardless whether the star is barotropic or non-barotropic. (Note that Lindblom & Ipser call any inertial mode 
a generalized r mode). 

Let us start with the barotropic case for m = 2. As already mentioned above, the axial equation for I = 2 yields the band of 
frequencies given by Eq. (|7^). This is depicted in the left graph of Fig. 3, where we show the continuous spectrum as a function 
of M/R. No inertial modes are present. We have also included the Newtonian value for the r mode as given by Eq. (^). Note 
that this value is completely independent of any stellar structure parameters. If we now include the coupling to the polar 
equations with I = 3, the situation changes drastically. We now obtain two continuous frequency bands and a mode which 
lies below the Newtonian r mode value, the further below the more compact the star is. Translating to positive frequencies, 
this means that this mode always has a higher frequency than the Newtonian r mode. Inspection of the eigenfunction reveals 
that the axial part fits nicely to a r l+1 power law, i.e. 

4 =2 ~ r\ (74) 

which is exactly what one expects for the r mode. Inclusion of higher I only gives minor frequency corrections, which tells us 
that this is indeed a physical inertial mode, the relativistic version of an m = 2 r mode. According to our nomenclature for 
inertial modes,the r mode is actually an i\ inertial mode. In Fig. 3, we have also included the coupling up to £ max = 4 and 
5. Here, each inclusion of one additional I leads to another patch in the continuous spectrum. Notice that the width of the 
uppermost and lowermost patches is very broad, whereas the middle patches around a /Ql ~ —1 remain very narrow, even 
for large M/R. For each £ max , we find £ max — m new inertial modes, appearing very close to those frequencies which were 
occupied by the continuous spectrum for £ max — 1. Again, some modes exist throughout the whole compactness range, while 
some are destroyed by the continuous spectrum above some M/R. In particular this is the case for those modes residing at 
both ends of the frequency range, as the continuous bands become very broad for large M/R. 

To assess the accuracy of our results, we have compared our r mode frequencies with non-linear time evolutions. For 
the 1.4 Mq neutron star model with fi — 2180 s" 1 , the non-linear evolution in the Cowling approximation yields an r-mode 
frequency of 500 Hz (Stergioulas, 2002). Our linear calculation gives a frequency of 512 Hz, which differs only by 2.5 per cent 
from the non-linear value. For larger rotation rates the difference increases, and for fi = 4986 s" 1 , the non-linear value is 
1030 Hz, whereas we find 1172 Hz, which is by 14 per cent off. 

When turning to the non-barotropic case the picture becomes more complicated because of the presence of the g modes. 
When the star is set into rotation, the g and inertial modes start to interact. Yoshida & Lee (2000b) have studied in depth the 
behaviour of the g and inertial modes as functions of the rotation rate and Ti for Newtonian stars. In this paper, we would 
like to focus only on the influence of the non-barotropicity on the r mode. 

In accordance with the results Yoshida & Lee (2000b) we find that the r mode frequency remains almost unaffected as 
the star starts to deviate from being barotropic. However, we find that for small rotation rates, the eigenfunctions change 
in a quite drastic way if Fi starts to differ from F. The larger the rotation, the more similar they get. We have already 
mentioned above that the relativistic r mode of a barotropic star has to be an axial-led hybrid mode, which must retain its 
polar contribution in the non-rotating limit, whereas for a non-barotropic star the polar contribution must vanish. This can be 
readily checked by computing a sequence of eigenfunctions for decreasing values of fi. These results are shown in Fig. 4, were 
we contrast the m = 2r modes for a barotropic model with those of the corresponding non-barotropic model with Ti = 3. 
We only show the dominant parts of the eigenfunctions, which are the axial contribution U3 with / = 2 and the polar parts 
H, Ui and U2 with / = 3. Furthermore, we rescale the eigenfunctions in such a way that u-z(R) = 1. 
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For the barotropic sequence (left panels), we can see that the shape of M3 does not depend on the rotation rate and is 
always given by r 3 . The polar parts H, u\ and 112 actually do change along the sequence. Whereas H goes to zero, ui and 
112 approach some finite values. This is in full agreement with Lockitch et al. (2001), who have shown that the non-trivial 
stationary perturbations must have vanishing H and non- vanishing ui and U2. From Eq. (^2|) is clear that for vanishing H 
(and therefore vanishing H), ui and 112 must satisfy 

' _I_ fo ,' \> J_ 2 \ 2x 1(1+1) 

Ui + I 2v — A + - j ui — e — -5 — u 2 = 0, (75) 

which is equivalent to Eq. (3.23) of Lockitch et al. (2001). 

The non-barotropic sequence (right panels) is completely different. For a large rotation rate e = 0.2, the eigenfunction is 
very similar to the barotropic one. However, for smaller rotation rates, they begin to deviate considerably. The shape of the 
axial part 113 is no longer described by r 3 , but now becomes zero throughout most parts of the star and only peaks close to 
the surface. Our chosen value of Ti = 3 seems to be quite large, however, the same effect can be observed for Fi much close 
to 2, in this case, it appear for much smaller rotation rates. In principle, we should not expect this behaviour, instead the 
shape of 113 should remain the same and only the polar perturbations should go to zero. From Fig. 4 we can see that latter is 
actually the case and H as well as ui and U2 all go to zero in the non-rotating limit. We can track this singular behaviour of 
M3 back to the frame dragging ui, for if we set it to zero, thus simulating a Newtonian stellar model, the eigenfunction of 113 
is well behaved; i.e. it remains proportional to r 3 for any rotation rate, whereas all the polar perturbations go to zero in the 
non-rotating limit. We suppose that this singular behaviour of 113 is an artefact of the Cowling approximation, which should 
be cured by taking into account the metric perturbations. 



4 SUMMARY 

We have investigated the inertial modes of slowly rotating relativistic stars in the Cowling approximation. Considering only 
first order rotational corrections, we have performed our computations in both the time and the frequency domain. Comparison 
with non-linear results for the /, p and r modes showed very good agreement, at least for slowly rotating stars. For larger 
rotation, second order effects become important and would have to be included for accurate mode calculations. 

For the numerical calculations, the system of equations has to be truncated at some cut-off value Z max . This truncation, 
however, introduces many spurious modes, which either vanish or shift in frequencies when £ max is increased. The "true" 
inertial modes, which appear for a given £ max are those which do not change significantly when i max becomes larger. For each 
/ > m > 0, one can find I — m inertial modes, which can be labeled i n (1 < n < I — m). For the axisymmetric case m = 0, 
there are / — 1 different modes, although in this case one mode has zero frequencies and the others always come in pairs with 
positive and negative frequencies. 

We have shown that, as a result of the relativistic frame dragging, there always exists a continuous spectrum, which, 
however, depends very strongly on the number of coupled equations. For Z max = 2, the continuous spectrum is confined to 
a single connected region, with an increasing width for more relativistic stellar models. If i max is increased the continuous 
spectrum splits into an increasing number of patches, which are disconnected for weakly relativistic stars but tend to overlap 
for more compact stellar models. For ra^O, each successive inclusion of a higher I creates a new patch, whereas for m = 0, 
only every second / does so. In the limit Z max — ■» 00 the continuous spectrum should completely cover the frequency range of 
the inertial modes. In previous work, which neglected the coupling between the equations, it was stated that modes should not 
be able to exist inside the continuous spectrum, as the eigenfunctions would become singular. Through the time evolutions, 
however, we were able to show that there are modes, which can exist inside the continuous spectrum. 

We provide the following (phenomenological) explanation. If an inertial mode can exist inside the continuous spectrum 
for a given i max , then it had to appear for the first time outside the continuous spectrum for some smaller i max . If a mode 
with label i l ° appears outside the continuous spectrum for Z max = lo, then it will not be affected if the continuous spectrum 
shifts itself on top of the mode for some £ max > lo. On the other hand, if the continuous spectrum for i max = lo is already so 
broad that it covers the position where the i l £ mode would appear, then this mode cannot come into existence. 

Hence, it seems that for very relativistic stars, the continuous spectrum is able to destroy some of the inertial modes. We 
should, however, be still very cautious to draw any conclusions about the (non)existence of these modes for a rapidly rotating 
neutron star. We have seen that the continuous spectrum is very sensitive to the number of equations that are coupled. As we 
are only taking into account the first order rotational corrections, there is only coupling from I to I ± 1. Inclusion of second 
order corrections leads to coupling to I ± 2; however, this should not affect the continuous spectrum, since in the relevant 
equations for U2 and 113, only terms coupling to H and ui arise. Third order corrections might do the job and it could well 
be that they will affect the continuous spectrum in such a way that modes which do not exist in the first order analysis can 
now exist because the responsible patch could have moved somewhere else. 

So far, we have completely neglected any metric perturbations, as we worked in the Cowling approximation. Lockitch 
et al. (2001) included some (non-radiative) metric perturbations, but restricted their studies to a post Newtonian treatment. 
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Nevertheless it is clear that we should expect some quantitative difference when including these terms, i.e. when we move from 
the Cowling approximation to the low-frequency approximation. Still, the continuous spectrum will remain unaffected when 
metric variables are included, for it is only the fluid equations which are responsible for its existence. By including the metric 
perturbations, the modes might be affected in such a way, however, that they get pushed out of the continuous spectrum and 
we might find some of the modes which we cannot find in the Cowling approximation. We think that this scenario seems to 
be quite likely as the study of the purely axial I — m case has shown. There, the Cowling approximation leads to a purely 
continuous spectrum without any mode solutions. In the low-frequency approximation, where one retains a certain metric 
component, one is lead to Kojima's master equation, which in certain cases admits r mode solutions, whose frequencies lie 
outside the continuous spectrum (Ruoff & Kokkotas 2001; Yoshida 2001). But the range of the continuous spectrum was still 
the same as in the Cowling approximation and even when all the metric perturbations were included (Ruoff & Kokkotas 
2002; Yoshida & Futamase 2001), the qualitative picture remained the same. Hence, we conjecture that the inclusion of the 
metric perturbations in the fully coupled equations will not affect the behaviour of the continuous spectrum, but it certainly 
will shift the modes, and possibly in such a way that some of them might be pushed outside the patches of the continuous 
spectrum and therefore reappear. 

As far as the r modes are concerned, our results show that their existence does not seem to called into question by 
the continuous spectrum. The previous results by Ruoff & Kokkotas (2001, 2002), Yoshida (2001) and Yoshida & Futamase 
(2002), which suggested that in certain stellar models with rather soft equations of state, the (purely axial) relativistic r 
modes cannot exist due to the existence of the continuous spectrum, seem to be artefacts arising from neglecting the coupling 
to the polar equations. The inclusion of this coupling shifts the location of the continuous spectrum and the r mode can 
exist well outside it for a large range of stellar models. In this paper, we have only presented results for a single poly tropic 
equation of state with F = 2, but from our previous results (Ruoff & Kokkotas 2001), we know that for V > 2, the width 
of the continuous spectrum actually shrinks and thus favours the presence of inertial modes. As most of the current realistic 
equations of state are stiffer than a F = 2 polytrope, we expect that the existence of r modes should not be questioned. The 
same should also hold for strange stars. 

Nevertheless it is clear that the relativistic case is quite different from the Newtonian one, and further studies are necessary 
to gain a full understanding of the inertial modes of (rapidly) rotating relativistic stars. In particular the metric perturbations 
have to be included to obtain the correct r-mode frequencies and the associated growth times. We cannot rule out that, in 
particular the latter will deviate significantly from the results obtained from applying the quadrupole formula to Newtonian 
stellar models. 
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Figure 1. Modes appearing for different cut-off values i max for a stellar model with M/R = 0.01. Around each patch of the continuous 
spectrum (grey shaded areas), there are many modes, which dissappear as Z max is increased. The "true" inertial modes, which are labeled 
with v n , remain fixed as i m ax is further increased. 
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Figure 3. Same as in Fig. 3 for m = 2. In addition we have plotted the Newtonian value for the r-mode according to Formula JF^). The 
relativistic r mode appears for i max = 3 well outside the continuous spectrum as the latter has shifted away from its previous location 
for i max = 2. Only the r mode and the and i| modes can exist for the whole compactness range. 
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Figure 4. Sequence of r-mode eigenfunctions for a barotropic (ri = V = 2, left panels) and a non-barotropic (Ti = 3, right panels) 
stellar model with M/R = 0.01. The plot style of the curves is solid for e = 0.2, long dashed for e = 0.1, short dashed for e = 0.05 and 
dotted for e = 0.01. 
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